Baseline dynamics
\[\begin{align}
S' & =
\mu N^\star - (\lambda + \mu) S
\\
E' & =
\lambda S - (\kappa E + \mu E)
\\
I_S' & =
p \kappa E - (\gamma_S + \delta_H + \mu_{I_S} + \mu) I_S
\\
I_A' &=
(1 - p) \kappa E - (\gamma_A + \mu) I_A
\\
H' &=
\delta_H I_S - (\gamma_H + \mu_H + \mu) H
\\
R' & =
\gamma_S I_S + \gamma_A I_A + \gamma_H H - \mu R
\\
D' &=
\mu_{I_S} I_S + \mu_H H
\\
\lambda &:=
\frac{\beta_A I_A + \beta_S I_S}{N^{\star}}
\\
N^{\star}(t) &=
S + E + I_S + I_A +
H + R
\end{align}\]
Observation Model
\[\begin{align}
Y_t &
\sim \text{Poisson}(\lambda_t), \qquad
\\
\lambda_t &
= \int_0^t \rho \delta_e E
\\
& p \sim \text{Uniform(0.3 0.8)}
\\
& \kappa: \sim \text{Gamma(10,50)}
\end{align}\]
Postulated densities
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
## Loading required package: StanHeaders
## rstan (Version 2.18.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
mu_ <- function(lambda, ylabel= "mu_h"){
# We suppose mu_. as a exponential r.v.
x <- seq(from = 0, to = 1, by=.001)
mu_h_ <- dexp(x, rate = lambda)
df <- data.frame(x=x, mu_h = mu_h_)
mean_ <- 1/lambda
mean_time <- 1/mean_
text <- paste(mean_time,"days")
p <- ggplot(data=df,
aes(x = x,
y = mu_h)) +
xlab("x")+
ylab(ylabel) +
geom_line(colour = "#000000") +
geom_area(alpha = 0.6, fill = "lightgray") +
geom_vline(xintercept = mean_, linetype="dotted",
color = "red", size=0.5)+
annotate(geom = "text", x = 1.30 * mean_, y=0, label = text,
color="black")
return(p)
}
#
gamma_ <- function(alpha, beta, ylabel="gamma_S"){
# We suppose mu_. as a exponential r.v.
theta <- 1 / beta
x <- seq(from = 0, to = 1, by=.001)
gamma_s_ <- dgamma(x, shape = alpha,
scale= theta)
df <- data.frame(x=x, gamma_s = gamma_s_)
mean_ <- alpha * theta
mean_time <- 1/mean_
text <- paste(mean_time,"days")
p <- ggplot(data=df,
aes(x = x,
y = gamma_s)) +
geom_line(colour = "#000000") +
geom_area(alpha = 0.6, fill = "lightgray") +
geom_vline(xintercept = mean_, linetype="dotted",
color = "red", size=0.5) +
annotate(geom="text", x=1.30 * mean_, y=0, label = text,
color="black")
return(p)
}
pi_mu_h <- mu_(lambda = 25, ylabel= "mu_h")
pi_mu_i_s <- mu_(lambda = 35, ylabel= "mu_i_s")
pi_gamma_s <- gamma_(alpha = 10, beta = 100, ylabel= "gamma_S")
pi_gamma_a <- gamma_(alpha = 10, beta = 50, ylabel= "gamma_A")
pi_gamma_h <- gamma_(alpha = 10, beta = 250, ylabel= "gamma_H")
fig6 <- ggplotly(pi_mu_h)
fig7 <- ggplotly(pi_mu_i_s)
fig8 <- ggplotly(pi_gamma_s)
fig9 <- ggplotly(pi_gamma_a)
fig10 <- ggplotly(pi_gamma_h)
fig6
Lockdown-Vaccination-Treatment-Model
\[\begin{align}
S'(t) &=
\mu \bar{N}
-\frac{\beta_S I_S+\beta_AI_A}{\bar{N}}S
-(\mu+\lambda_V)S +\delta_V V
\\
E'(t) &=
\frac{\beta_S I_S
+ \beta_A I_A}{\bar{N}} S
+ \epsilon \frac{\beta_S I_S + \beta_A I_A}{\bar{N}}V
- (\mu+\delta_E) E \nonumber
\\
I'_S(t) & =
p \delta_E E
-(\mu + \mu_S + \alpha_S + \lambda_T) I_S
+(1 - q) \alpha_T T
\\
I'_A(t) & =
(1-p) \delta_E E
- (\mu + \mu_A + \alpha_A) I_A
\\
R'(t) & =
\alpha_S I_S
+ \alpha_A I_A
+ q \alpha_T T
- \mu R
\\
D'(t)& =
\mu_S I_S
+ \mu_A I_A
\\
V'(t)& =
\lambda_V S
- \epsilon
\frac{\beta_S I_S + \beta_A I_A}{\bar{N}} V
- (\mu + \delta_V) V
\\
T'(t) &=
\lambda_T I_S
-(\mu + \alpha_T) T
\\
\bar{N}(t) &=
S(t) + E(t) + I_S(t) + I_A(t) + R(t) + V(t) + T(t)
\end{align}\]
Stan MODEL
## functions {
## real[] seirvt(real t, // time
## real[] y,
## real[] theta,
## real[] x_r,
## int[] x_i) {
## //
## real dy_dt[8];
## real s = y[1];
## real e = y[2];
## real i_s = y[3];
## real i_a = y[4];
## real h = y[5];
## real r = y[6];
## real d = y[7];
## real n_star;
## //
## real beta_s = theta[1];
## real beta_a = theta[2];
## real kappa = theta[3];
## real p = theta[4];
## real delta_h = theta[5];
## real mu_i_s = theta[6];
## real mu_h = theta[7];
## real gamma_s = theta[8];
## real gamma_a = theta[9];
## real gamma_h = theta[10];
## //
## real force_infection;
## real mu = 3.913894e-05;
## n_star = s + e + i_s + i_a + h + r;
##
## force_infection = (beta_s * i_s + beta_a * i_a) / n_star;
## dy_dt[1] = mu * n_star - force_infection * s - mu * s;
## dy_dt[2] = force_infection * s - (mu + kappa) * e;
## dy_dt[3] = p * kappa * e - (gamma_s + mu_i_s + delta_h) * i_s;
## dy_dt[4] = (1.0 - p) * kappa * e - (gamma_a + mu) * i_a;
## dy_dt[5] = delta_h * i_s - (gamma_h + mu_h + mu) * h;
## dy_dt[6] = gamma_s * i_s + gamma_a * i_a + gamma_h * h - mu * r;
## dy_dt[7] = mu_i_s * i_s + mu_h * h;
## dy_dt[8] = p * kappa * i_s;
## return dy_dt;
## }
## }
## data {
## int<lower = 1> n_obs; // number of days observed
## int<lower = 1> n_theta; // number of model parameters
## int<lower = 1> n_difeq; // number of differential equations
## int<lower = 1> n_pop; // population
## int y[n_obs]; // data, total number of infected
## real t0; // initial time point (zero)
## real ts[n_obs]; // time points observed
## }
##
## transformed data {
## real x_r[0];
## int x_i[0];
## }
## parameters {
## real <lower = 5e-6> theta[n_theta];
## // real <lower = 0, upper = 1> E0;
## // real <lower = 0, upper=1> IA0;
## real <lower = 0, upper = n_pop> E0;
## real <lower = 0, upper = n_pop> IA0;
## }
## transformed parameters{
## real y_hat[n_obs, n_difeq]; // solution from the ODE solver
## real y_init[n_difeq];
## // initial conditions
## // real IS0 = 8.297217e-06;
## real IS0 = 22.0;
## real H0 = 0.0;
## real R0 = 0.0;
## real D0 = 0.0;
## real CIS0 = 74.0;
## y_init[1] = n_pop - (IS0 + IA0 + E0);
## y_init[2] = E0;
## y_init[3] = IS0;
## y_init[4] = IA0;
## y_init[5] = H0;
## y_init[6] = R0;
## y_init[7] = D0;
## y_init[8] = CIS0;
## y_hat = integrate_ode_rk45(seirvt, y_init, t0, ts, theta, x_r, x_i);
## }
## model {
## real lambda[n_obs]; // poisson parameter
## // priors
## theta[1] ~ normal(1.0, 0.1); // beta_s
## theta[2] ~ normal(1.0, 0.1); // beta_a
## theta[3] ~ gamma(10, 40); // kappa
## theta[4] ~ uniform(0.1, 0.5); // p
## theta[5] ~ gamma(10, 40); // delta_h
## theta[6] ~ exponential(33.33333); // mu_i_s
## theta[7] ~ exponential(25); // mu_h
## theta[8] ~ gamma(10, 100); // gamma_s
## theta[9] ~ gamma(10, 50); // gamma_a
## theta[10] ~ gamma(10, 250); // gamma_h
## //
## // E0 ~ uniform(1.121246e-07, 3.363737e-05);
## // IA0 ~ uniform(7.848719e-06, 3.363737e-05);
## E0 ~ uniform(74, 2100);
## IA0 ~ uniform(74, 2100);
## // likelihood
## for (i in 1:n_obs){
## //lambda[i] = y_hat[i, 8] * n_pop;
## lambda[i] = y_hat[i, 8];
## }
## y ~ poisson(lambda);
## }
## generated quantities {
## real R_0; // Basic reproduction number
## real mu = 3.913894e-05;
## real beta_s = theta[1];
## real beta_a = theta[2];
## real p = theta[4];
## real gamma_s = theta[4];
## real gamma_a = theta[5];
## real kappa = theta[6];
## R_0 = p * beta_s * kappa / ((gamma_s + mu) * (kappa + mu))
## + (1.0 - p) * beta_a * kappa / ((gamma_a + mu) * (kappa + mu));
## }
Observation Model
\[\begin{align}
Y_t &
\sim \text{Poisson}(\lambda_t), \qquad
\\
\lambda_t &
= \int_0^t p \kappa E
\\
& p \sim \text{Uniform(0.3 0.8)}
\\
& \kappa \sim \text{Gamma(10, 52)}
\end{align}\]